The inner product, also known as the dot product can be computed in R using the following function:
Given the two vectors:
u <- c(2, -5, -1)
v <- c(3, 2, -3)
t(u) %*% v
## [,1]
## [1,] -1
dot_product <- function(vector1, vector2) {
t(vector1) %*% vector2
}
dot_product(u, v)
## [,1]
## [1,] -1
The length (or norm) of a vector is the square root the sum of all entries in the vector squared.
A vector whose length is 1 is called a unit vector. Example:
Let v = (1, -2, 2, 0). Find a unit vector u in the same direction as v.
Compute the length of v:
v <- c(1, -2, 2, 0)
length_vector <- function(vector, norm = FALSE) {
length_v <- sqrt(sum(vector*vector))
## Multiply __v__ by 1/length_v
if (norm){
return((1/length_v)*vector)
}
else {
return(length_v)
}
}
#Length of a vector:
length_vector(v)
## [1] 3
#Unit Vector
length_vector(v, TRUE)
## [1] 0.3333333 -0.6666667 0.6666667 0.0000000
# Check that vector is a unit vector
# It returns 1 so this is true.
length_vector(length_vector(v, TRUE))
## [1] 1
For u and v in Rn, the distance between u and v, written as dist(u, v), is the length of the vector u-v. That is, dist(u, v) = ||u - v||
Example: Compue the distance between the vectors u = (7,1) and v = (3,2)
# Let
u <- c(7, 1)
v <- c(3, 2)
u_minus_v <- u - v
length_vector(u_minus_v)
## [1] 4.123106
Two vectors are orthogonal to each other if the dot product is 0.
u <- c(0, 1)
v <- c(1, 0)
dot_product(u, v)
## [,1]
## [1,] 0
Orthogonal Sets
u1 <- c(3, 1, 1)
u2 <- c(-1, 2, 1)
u3 <- c(-1/2, -2, 7/2)
dot_product(u1, u2)
## [,1]
## [1,] 0
dot_product(u2, u3)
## [,1]
## [1,] 0
dot_product(u1, u3)
## [,1]
## [1,] 0
An orthogonal basis for a subsace W of Rn is a basis for W that is also an orthogonal set.
y <- c(6, 1, -8)
# Returns y
dot_product(y, u1)/dot_product(u1, u1)*u1+dot_product(y, u2)/dot_product(u2, u2)*u2 + dot_product(y, u3)/dot_product(u3, u3)*u3
## [1] 6 1 -8
Orthogonal projections
y <- c(7, 6)
u <- c(4, 2)
orthogonal_proj <- function(y, Un, y_minus_yhat = FALSE) {
y_hat <- dot_product(y,Un)/dot_product(Un,Un)*Un
if(y_minus_yhat) {
return(y - y_hat)
}
else {
return(y_hat)
}
}
# Orthogonal projection of y onto u
orthogonal_proj(y, u)
## [1] 8 4
# Component of y orthogonal to u
orthogonal_proj(y, u, TRUE)
## [1] -1 2
# Verify that they sum to y
orthogonal_proj(y, u) + orthogonal_proj(y, u, TRUE)
## [1] 7 6
# Verify that they are perpindicular
dot_product(orthogonal_proj(y, u), orthogonal_proj(y, u, TRUE))
## [,1]
## [1,] 0
length_vector(orthogonal_proj(y, u, TRUE))
## [1] 2.236068
v1 <- c(3/sqrt(11), 1/sqrt(11), 1/sqrt(11))
v2 <- c(-1/sqrt(6), 2/sqrt(6), 1/sqrt(6))
v3 <- c(-1/sqrt(66), -4/sqrt(66), 7/sqrt(66))
dot_product(v1, v2)
## [,1]
## [1,] 0
dot_product(v1, v3)
## [,1]
## [1,] 5.551115e-17
dot_product(v2, v3)
## [,1]
## [1,] 0
dot_product(v1, v1)
## [,1]
## [1,] 1
dot_product(v2, v2)
## [,1]
## [1,] 1
dot_product(v3, v3)
## [,1]
## [1,] 1
U <- matrix(c(v1, v2, v3), ncol = 3)
round(t(U)%*%U, 1)
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
# I used round to show that this is nearly the identity matrix.
t(U)%*%U
## [,1] [,2] [,3]
## [1,] 1.000000e+00 7.681559e-19 3.007996e-17
## [2,] 7.681559e-19 1.000000e+00 2.696773e-17
## [3,] 3.007996e-17 2.696773e-17 1.000000e+00
u1 <- c(2, 5, -1)
u2 <- c(-2, 1, 1)
y <- c(1, 2, 3)
# The orthogonal projection of y onto W
orthogonal_proj(y, u1) + orthogonal_proj(y, u2)
## [1] -0.4 2.0 0.2
y - (orthogonal_proj(y, u1) + orthogonal_proj(y, u2))
## [1] 1.4 0.0 2.8
u1 <- c(2, 5, -1)
u2 <- c(-2, 1, 1)
y <- c(1, 2, 3)
# W = Span {u1, u2}
# The closest point in W to y is
orthogonal_proj(y, u1) + orthogonal_proj(y, u2)
## [1] -0.4 2.0 0.2
y <- c(-1, -5, 10)
u1 <- c(5, -2, 1)
u2 <- c(1, 2, -1)
# The distance from y to W is ||y-y_hat||
y - orthogonal_proj(y, u1) - orthogonal_proj(y, u2)
## [1] 0 3 6
length_vector(y - orthogonal_proj(y, u1) - orthogonal_proj(y, u2))
## [1] 6.708204
[](
library(pracma)
x1 <- c(1, 1, 1, 1)
x2 <- c(0, 1, 1, 1)
x3 <- c(0, 0, 1, 1)
# Construct an orthogonal basis for W
# Step 1:
v1 <- x1
# Step 2:
v2 <- x2 - orthogonal_proj(x2, v1)
v1
## [1] 1 1 1 1
v2
## [1] -0.75 0.25 0.25 0.25
#Step 3
v3 <- x3 - orthogonal_proj(x3, v1) - orthogonal_proj(x3, v2)
v3
## [1] 0.0000000 -0.6666667 0.3333333 0.3333333
A <- matrix(c(x1, x2, x3), ncol = 3)
gS <- gramSchmidt(A)
round(gS$Q %*% gS$R, 1)
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 1 1 0
## [3,] 1 1 1
## [4,] 1 1 1
v1 <- c(3, 6, 0)
v2 <- c(0, 0, 2)
u1 <- v1/length_vector(v1)
u2 <- v2/length_vector(v2)
u1
## [1] 0.4472136 0.8944272 0.0000000
u2
## [1] 0 0 1
v1 <- c(1, 1, 1, 1)
v2 <- c(-3, 1, 1, 1)
v3 <- c(0, -2/3, 1/3, 1/3)
u1 <- v1/length_vector(v1)
u2 <- v2/length_vector(v2)
u3 <- v3/length_vector(v2)
Q <- matrix(c(u1, u2, u3), ncol = 3)
Q
## [,1] [,2] [,3]
## [1,] 0.5 -0.8660254 0.00000000
## [2,] 0.5 0.2886751 -0.19245009
## [3,] 0.5 0.2886751 0.09622504
## [4,] 0.5 0.2886751 0.09622504
A
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 1 1 0
## [3,] 1 1 1
## [4,] 1 1 1
R <- t(Q) %*% A
R
## [,1] [,2] [,3]
## [1,] 2.000000e+00 1.5000000 1.0000000
## [2,] 1.110223e-16 0.8660254 0.5773503
## [3,] 0.000000e+00 0.0000000 0.1924501
round(R, 7)
## [,1] [,2] [,3]
## [1,] 2 1.5000000 1.0000000
## [2,] 0 0.8660254 0.5773503
## [3,] 0 0.0000000 0.1924501
gramSchmidt(A)
## $Q
## [,1] [,2] [,3]
## [1,] 0.5 -0.8660254 0.0000000
## [2,] 0.5 0.2886751 -0.8164966
## [3,] 0.5 0.2886751 0.4082483
## [4,] 0.5 0.2886751 0.4082483
##
## $R
## [,1] [,2] [,3]
## [1,] 2 1.5000000 1.0000000
## [2,] 0 0.8660254 0.5773503
## [3,] 0 0.0000000 0.8164966
A <- matrix(c(4, 0, 0, 2, 1, 1), ncol = 2, byrow = TRUE)
b <- c(2, 0, 11)
solve(t(A)%*%A)%*%t(A)%*%b
## [,1]
## [1,] 1
## [2,] 2
least_squares_sol <- function(A, b) {
solve(t(A)%*%A)%*%t(A)%*%b
}
least_squares_sol(A, b)
## [,1]
## [1,] 1
## [2,] 2
# Least squares error:
b <- c(2, 0, 11)
A <- matrix(c(4, 0, 0, 2, 1, 1), ncol = 2, byrow = TRUE)
x_hat <- c(1, 2)
length_vector(b - A%*%x_hat)
## [1] 9.165151
A <- matrix(c(1, 3, 5, 1, 1, 0, 1, 1, 2, 1, 3, 3), ncol = 3, byrow = TRUE)
b <- c(3, 5, 7, -3)
least_squares_sol_gs <- function(A, b) {
gS <- gramSchmidt(A)
solve(gS$R)%*%t(gS$Q)%*%b
}
least_squares_sol_gs(A, b)
## [,1]
## [1,] 10
## [2,] -6
## [3,] 2
least_squares_sol(A, b)
## [,1]
## [1,] 10
## [2,] -6
## [3,] 2
# Fit a line
# Two arbitrary vectors of equal length
x <- c(0, 1, 2, 3, 4, 5, -2, 3, 5)
y <- c(0, 1, 4, 9, 16, 25, 0, 3, 7)
# Matrix Computation
A <- matrix(c(rep(1, length(x)),x), ncol = 2)
A
## [,1] [,2]
## [1,] 1 0
## [2,] 1 1
## [3,] 1 2
## [4,] 1 3
## [5,] 1 4
## [6,] 1 5
## [7,] 1 -2
## [8,] 1 3
## [9,] 1 5
AtA <- t(A)%*%A
B <- t(A) %*% y
vector <- solve(AtA)%*%B
vector
## [,1]
## [1,] 1.000000
## [2,] 2.666667
# Plot points and line
slope <- vector[2,]
intercept <- vector[1,]
plot(x, y)
curve(intercept+slope*x, add = TRUE)
# Fit a parabola
# Two arbitrary vectors of equal length
x <- c(0, 1, 2, 3, 4, 5, -2, 3, 5)
y <- c(0, 1, 4, 9, 16, 25, 0, 3, 7)
# Matrix Computation
A <- matrix(c(rep(1, length(x)),x, x^2), ncol = 3)
A
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 1 1 1
## [3,] 1 2 4
## [4,] 1 3 9
## [5,] 1 4 16
## [6,] 1 5 25
## [7,] 1 -2 4
## [8,] 1 3 9
## [9,] 1 5 25
AtA <- t(A)%*%A
B <- t(A) %*% y
vector <- solve(AtA)%*%B
vector
## [,1]
## [1,] -0.0998308
## [2,] 0.9949239
## [3,] 0.4839255
# Plot points and line
slope <- vector[2,]
intercept <- vector[1,]
squared <- vector[3,]
plot(x, y)
curve(intercept+slope*x+squared*x*x, add = TRUE)
# Fit anything?
# Two arbitrary vectors of equal length
x <- c(0, 1, 2, 3, 4, 5, -2, 3, 5)
y <- c(0, 1, 4, 9, 16, 25, 0, 3, 7)
# Matrix Computation
A <- matrix(c(rep(1, length(x)),x, x^2), ncol = 3)
A
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 1 1 1
## [3,] 1 2 4
## [4,] 1 3 9
## [5,] 1 4 16
## [6,] 1 5 25
## [7,] 1 -2 4
## [8,] 1 3 9
## [9,] 1 5 25
AtA <- t(A)%*%A
B <- t(A) %*% y
vector <- solve(AtA)%*%B
vector
## [,1]
## [1,] -0.0998308
## [2,] 0.9949239
## [3,] 0.4839255
# Plot points and line
slope <- vector[2,]
intercept <- vector[1,]
squared <- vector[3,]
plot(x, y)
curve(intercept+slope*x+squared*x*x, add = TRUE)
X <- matrix(c(1, 2, 1, 5, 1, 7, 1, 8), ncol = 2, byrow = TRUE)
y <- c(1, 2, 3, 3)
least_squares_line <- function(X, y) {
coefficients <- solve(t(X)%*%X)%*%t(X)%*%y
plot(X[,2], y)
curve(coefficients[1,] + coefficients[2,]*x, add = TRUE)
}
least_squares_line(X, y)
#Playing with some dui data
dui_month <- read.csv("C:/Users/Freddy/Documents/dui_month.csv", stringsAsFactors = FALSE, header = TRUE)
dui_month <- dui_month[,-1]
dui_month <- t(dui_month)
dui_month
## [,1] [,2] [,3] [,4] [,5]
## Jan 464 459 548 502 555
## Feb 421 431 559 514 509
## Mar 455 573 609 592 484
## Apr 412 388 519 522 481
## May 418 495 539 553 445
## Jun 352 527 533 484 474
## Jul 395 489 512 486 479
## Aug 393 547 547 539 542
## Sep 449 500 491 441 484
## Oct 451 511 499 539 525
## Nov 368 527 488 492 481
## Dec 487 560 525 491 385
length(dui_month)
## [1] 60
plot(1:60, dui_month)
dui_month
## [,1] [,2] [,3] [,4] [,5]
## Jan 464 459 548 502 555
## Feb 421 431 559 514 509
## Mar 455 573 609 592 484
## Apr 412 388 519 522 481
## May 418 495 539 553 445
## Jun 352 527 533 484 474
## Jul 395 489 512 486 479
## Aug 393 547 547 539 542
## Sep 449 500 491 441 484
## Oct 451 511 499 539 525
## Nov 368 527 488 492 481
## Dec 487 560 525 491 385
dui <- matrix(dui_month, ncol = 1, byrow = FALSE)
plot(1:60, dui)
X <- matrix(c(rep(1, 60), 1:60, (1:60)^2, (1:60)^3), ncol = 4)
least_squares_line <- function(X, y) {
coefficients <- solve(t(X)%*%X)%*%t(X)%*%y
plot(X[,2], y)
curve(coefficients[1,] + coefficients[2,]*x + coefficients[3,]*x^2 + coefficients[4,]*x^3, add = TRUE)
}
least_squares_line(X, dui)